from pyscf import gto,lib,scf
import numpy as np
from opt import cre_coord

def cal_hess(coord,natom,basis,atom):
    hess=np.empty((3*natom,3*natom))
    count=0 
    for i in range(natom):
        for j in range(3):
            coord_m1=coord.copy()
            coord_m1[i,j] -= 1e-4 
            coord_p1=coord.copy()
            coord_p1[i,j] += 1e-4
            grad_m1 =scf.RHF(cre_coord(coord_m1,atom,basis)).run().nuc_grad_method().run().de
            grad_p1 =scf.RHF(cre_coord(coord_p1,atom,basis)).run().nuc_grad_method().run().de
            grad_tem = (grad_p1 - grad_m1) / 2e-4  ##grad_tem array:(natom,3)
            hess[i*3+j,:]=grad_tem.flatten()
            count +=1
    print('ok')
    return hess
